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Abstract 


A time  dependent,  three-dimensional  numerical  dynamic  model  of  the 
North  Pacific  Ocean,  possessing  an  actual  coastal  configuration  and  ten 


layers,  has  been  developed  in  order  to  show  the  physical  nature  of  the 
large-scale  normal  and  abnormal  characteristics  of  the  ocean  in  response 
to  the  various  normal  and  anomalous  seasonal  meteorological  conditions. 
Based  on  the  simulated  energetics,  emphasis  is  given  to  identify  the  major 
physical  processes  and  essential  dynamic  mechanisms  responsible  for  the 
generation,  evolution,  and  dissipation  of  large-scale  anomalies  in  the 
North  Pacific  Ocean.  The  model  is  based  on  time  integrations  of  the  fi- 
nite difference  forms  of  the  primitive  equations.  The  oceanic  circula- 
tion is  driven  by  atmospheric  forcing,  namely,  the  surface  wind  stresses 
and  the  differential  heating  over  the  ocean.  The  flux  form  of  numerical 
scheme  for  energy  conservation  and  the  rigid-lid  approximation  for  filter- 
ing out  the  external  inertia-gravity  wave  are  used  in  the  formulation. 

The  model  has  run  for  more  than  80  years,  simulated  with  the  annual 
mean  atmospheric  data  as  the  forcing  boundary  conditions.  The  long-term 
mean  state  in  the  model  reveals  the  large-scale  features  of  the  circulation 

patterns  and  density  distributions  in  the  North  Pacific  Ocean.  Three 

K 

\ N . 

gyres,  one  large  anticyclonic  in  the  subtropical  region  and  two  smaller 
ones  in  the  subartic  and  in  the  tropic  region,  are  well  developed.  The 
total  transport  near  the  western  boundary  reached  56  Sverdrups,  which 
agrees  quantitatively  well  with  observed  mean  transport  in  the  Kuroshio 
Current  south  of  Japan.  All  the  equatorial  currents,  the  Oyashio  Current. 
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the  Alaska  Current,  the  west  wind  drift,  and  the  California  Current  are 
grossly  shown.  Upwelling  is  observed  along  both  the  equator  and  the 
western  boundary.  Temperature  and  salinity  distributions,  except  in  high 
latitudes,  compare  well  with  observational  atlas  and  field  data. 

Further  experiments  for  the  time-dependent,  season  fluctuations  and 
for  the  thermal  anomaly  dynamics  in  the  North  Pacific  Ocean  have  been 
carried  out  and  will  be  reported  in  succeeding  papers. 


1.  Introduction 


In  recent  years,  the  substantial  deviations  of  sea-surface  tempera- 
tures from  their  long-term  means,  i.e.,  the  thermal  anomalies  in  the 
ocean,  have  attracted  the  attention  of  many  oceanographers  and  meterolo- 
gists.  These  anomalous  phenomena,  together  with  their  temporal  correla- 
tions and  their  spatial  teleconnections  with  oceanic  and  atmospheric 
changes,  have  been  extensively  studied  by  Bjerknes  (1962,  196 6,  1969, 

1972),  by  Namias  (1959,  1965,  1969,  1970),  and  by  others  (e.g.,  Petterssen, 
et  al. , 1962;  Rodewald,  1963).  The  general  conclusions  ascertain  that 
the  thermal  anomalies  in  the  ocean  are  constantly  modifying  the  overlying 
atmospheric  circulation,  thereby  altering  the  wind  systems  and  the  weather 
patterns.  Observations  indicate  that  climatic  fluctuations  in  the  ocean 
and  atmosphere  occur  in  all  parts  of  the  world  at  all  times  (Namias,  1965). 
In  the  Pacific  Ocean,  field  data  from  the  North  Pacific  Experiment  (NORPAX) 
have  shown  that  large  temperature  anomalies  with  intensities  sometimes  as 
high  as  2°C  to  5°C  frequently  occupy  one-third  to  one-half  or  more  of  the 
North  Pacific  Ocean,  and  they  usually  persist  for  long  periods  of  time. 
Besides  their  possible  feedback  to  the  atmosphere,  the  anomalies  cause  ab- 
normal advections  which,  together  with  the  atmospheric  anaomalous  varia- 
tions, influence  the  general  oceanic  circulation  patterns  (Namias,  1965, 
Arthur,  1966).  For  example,  in  1957-1958  an  anomalous  sea  surface  tempera- 
ture in  the  Pacific  was  associated  with  the  weakening  of  the  California 
Current,  and  forced  it  to  swing  offshore  (Isaacs,  1968).  This  resulted 
in  a strong  invation  of  warm,  subtropical  water  into  the  west  coast  as  far 


k 


r 
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north  as  British  Columbia. 

The  pliysicnl  natures  of  the  atmosphere,  the  ocean,  and  the  biota 
therein  are  closely  coupled  and  mutually  interacting.  In  the  large- 
scale  air/sea  interactions  the  ocean  primarily  supplies  heat  to  the  atmo- 
sphere while  the  atmosphere  supplies  momentum  and  heat  to  the  ocean.  The 
heat  released  from  the  ocean  modifies  the  atmospheric  circulation  while 
the  momentum  and  heat  exchange  from  the  atmosphere  to  the  ocean  determine 
the  flow  pattern  and  heat  content  in  the  ocean.  The  oceanic  circulation, 
then,  redistributes  the  thermal  energy,  thereby  setting  up  a closed  feed- 
back system  between  the  ocean  and  the  atmosphere.  For  the  purpose  of 
studying  the  mutual  effects  and  their  teleconnections  in  the  ocean/ 
atmosphere  system,  it  is  indispensible  to  employ  an  ocean/atmosphere 
coupled  system.  In  a pioneering  work,  Manabe  (1969)  and  Bryan  (1969a) 
performed  the  first  numerical  experiment  of  an  ocean/atmosphere  model 
that  included  major  coupling  and  mutual  adjusting  mechanisms.  Their  re- 
sults showed  that  the  poleward  transport  of  heat  by  oceans  reduced  the 
atmosphereic  meridional  heat  transport,  which  in  turn  decreased  the  baro- 
clinic  instability  in  the  troposphere.  In  the  tropics,  the  oceanic  ef- 
fect changed  the  rainfall  distribution  and  influenced  the  strength  of 
the  Hadley  circulation  producing  more  realistic  results.  Therefore,  the 
coupled  model  did  demonstrate  its  ability  to  more  accurately  portray  the 
large-scale  features  in  the  atmosphere. 

Simulation  studies  of  the  ocean  have  been  successfully  carried  out 
by  Takano  (1969),  II' in,  et  al.  (1969),  O’Brien  (1971),  Bye  and  Sag 
(1972),  Bryan  and  Cox  (1972),  and  Holland  and  Hirschman  (1975).  These 


5 


have  taken  into  account  realistic  configurations  of  the  basins  under 
their  investigation  with  a simplified  formulation  of  the  systems.  Bryan 
(1963,  1969b),  Gormatyuk  and  Sarkisyan  (1965),  Bryan  and  Cox  (1967,  1968), 
Crowley  (1970),  and  Haney  ( 197*0  have  taken  the  complete  dynamics  into 
their  models  of  certain  idealized  domains.  Cox  (1970)  has  used  Bryan's 
model  (1969b),  which  includes  the  whole  dynamics  as  well  as  the  irregu- 
lar boundaries,  to  investigate  the  monsoon  effect  in  the  circulation 
patterns  in  the  Indian  Ocean. 

Though  it  is  not  realistic  to  employ  a regional  atmospheric  model 
for  studies  of  planetary  responses  and  global  teleconnections,  a regional 
ocean  model  is  tolerable  for  dynamic  studies  in  the  ocean  due  to  its 
solid  boundary  enclosement  and  the  semi-isolated  characteristics  of  cer- 
tain basins.  If  one  addresses  the  more  specific  problems  of  physical 
processes,  thermal  structures  and  dynamic  responses  in  an  ocean  under 
meteorological  forcing,  it  is  probably  feasible  to  examine  a semi- 
enclosed  domain  such  as  the  North  Pacific  Ocean. 

In  the  investigation  of  the  anomalous  state  in  the  ocean,  the  an- 
omalies can  be  considered  as  the  perturbation  quantities,  though  not 
necessarily  very  small,  superimposed  upon  the  basic  long-term  averages. 

In  general,  the  oceanic  anomalies  are,  at  least  to  the  lower  order  of 
approximation,  generated  by  the  atmospheric  anomalies,  such  as  the  an- 
omalous pressure,  the  anomalous  wind,  the  anomalous  heat-exchange  and  water 
flux.  An  approximate  method  for  estimating  the  large-scale  sea-surface 
temperature  anomalies  due  to  atmospheric  changes,  proposed  by  Namias 
(1965)  and  extended  by  Arthur  (1 966)  and  Jacobs  (1967),  has  shown  some 


degree  of  success. 

We  have  investigated  the  anomaly  dynamics  in  the  North  Pacific  Ocean 
by  numerical  simulation.  The  ocean  model  we  developed  for  the  anomaly 
study  considers  the  actual  configuration  of  the  North  Pacific  Ocean  and 
takes  the  complete  dynamic  equations  with  all  large-scale  physical  pro- 
cesses into  consideration.  The  basic  numerical  scheme  used  in  the  model 
is  a modification  of  Haney's  idealized  ocean  model  (197*0 . The  space- 
staggered  grid  and  the  spacing  differencing  scheme  used  in  the  present 
model  follows  that  of  the  UCLA,  atmospheric  general  circulation  model 
(Gates,  et  ad.,  1971).  In  addition,  the  present  model  treats  the  upper 
boundary  conditions  in  an  improved  manner.  The  variation  of  heat  flux 
between  the  ocean  and  the  atmosphere  is  not  only  a function  of  the 
coupled  air-sea  temperatures  but  also  a summation  of  all  heating  and 
interacting  components  calculated  from  empirical  formulae  for  the  heat 
balance  at  the  naviface.  In  computing  wind  stress,  the  near  surface 
atmospheric  stability  is  also  taken  into  consideration.  Both  the 
longitudinal  component  and  the  meridional  component  of  wind  stress  are 
used  as  the  boundary  conditions  in  the  momentum  equation.  Furthermore, 
since  it  is  well-known  that  salinity  plays  an  important  role  in  the 
ocean  dynamics  at  high  latitudes,  a prognostic  equation  for  salinity 
conservation  is  also  included  in  the  hydrodynamic  system. 

Our  study  is  originated  from  the  objective  of  NORPAX  to  understand 
the  low-frequency  fluctuations  in  the  upper  layer  of  the  North  Pacific 
Ocean  and  their  relations  to  the  overlying  atmosphere.  The  major  goal 
is  to  investigate  the  anomaly  dynamics  and  to  describe  the  configuration 
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and  evolution  of  the  thermal  anomalies  in  the  North  Pacific  by  numerical 
simulations.  Because  the  anomalies  are  generally  small,  it  is  necessary 
to  take  a systematic  step-by-step  approach,  from  the  spin-up  of  an  ocean 
model  to  the  final  anomaly  experiments.  A set  of  three  papers  are  in 
progress  to  report  the  sequence  of  numerical  dynamic  study  for  the  North 
Pacific  Ocean.  This  first  paper  describes  briefly  the  basic  physics  in 
the  numerical  model  which  is  the  basic  tool  for  our  thermal  anomaly  in- 
vestigations, and  shows  the  simulated  long-term  mean  state  in  the  North 
Pacific  Ocean.  Our  major  interest  is  in  the  large-scale  (more  than 
1000  km)  and  long-term  (months  to  years)  fluctuations  in  the  North  Pacific 
Ocean  in  relation  to  the  overlying  atmosphere.  The  following  paper,  II, 
will  present  the  seasonal  fluctuations  in  the  North  Pacific  Ocean,  and 
III  will  discuss  the  dynamic  mechanisms  and  results  of  anomaly  experi- 
ments based  on  real-time  specifications  in  the  North  Pacific  model. 


* 
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2.  The  dynamic  model 


The  formulation  of  the  model  is  based  on  the  primitive  equations 
relative  to  the  right-hand  spherical  coordinate  system  (\,$,z)  on  the 
surface  of  the  earth  with  z taken  positive  upward  from  the  mean  sea  level. 
The  Navier-Stokes  equations  for  incompressible  and  nonhomogeneous  fluid 
with  some  basic  physical  approximations,  namely,  the  Boussinesq  approxi- 
mation, the  hydrostatic  approximation,  the  eddy  viscosity  approximation, 
and  the  traditional  approximation  (Eckart,  i960),  are  integrated  with 
time.  Let  (u,v)  be  the  horizontal  component  of  velocity  \V  and  define 
an  advective  flux  operator  as : 


(q)  = [(Uq)  + (vqcos4»)  ] + (wq) 


where  q is  any  scalar  quantity.  With  all  symbols  as  defined  in  Table  1, 
the  governing  equations  are: 


u = ^ u + vsin#(2fi+a  usec<>)  - (p  a)  ^sec^p  + X 
■fc  o \ 


v^  = v - usin<t>(2ft+a  ^usecJ>)  - (p^a)  + X 


pz  = - pg 


a sec<>[u^+(  vcos<t)  + w^  = 0 


T = >fT  + Q 

X 


s = >r  s + h 

X 


p = p(T,S,p) 


where  all  alphabetic  subscripts  denote  partial  derivatives  and 


II 

2-2  2 -2 

v,  [V  u+a  ( 1-tan  <t> ) u - 2a  tan'&sin'l’v  ] + v n 

1 \ 2 zz 

(if) 

$ 

F 

2-2  2 -2 

v^[V  v+a  ( 1-tan  <t>)v  + 2a  tan<Dsin<t>u^]  + VgV 

(10) 

2 -1 

(11) 

Q = 

k. V T + h klT 

1 2 zz 

2 -1 

(12) 

H = 

k V S + h k S 

1 2 zz 

wherein 


2 -2  2 -1  -k 

V (q)  = a sec  4>q^  + a sec<l>(a  cos<t>q^)( 


(13) 


Following  Bryan  (1969b),  the  step  function  h(-p^)  in  (ll)  and  (12)  incor- 
porates the  convective  density  adjust  mechanism  wherein  p ' indicates  the 
potential  density  of  the  ocean  water.  Notice  that,  although  in  most 
cases  the  vertical  column  of  the  oceanic  water  is  stable,  static  insta- 
bility due  to  inversions  in  the  potential  density  stratification  can  occur 
at  any  level  and  cause  strong  vertical  mixing  in  the  ocean.  Nevertheless, 
most  of  the  unstable  phenomena  happen  in  the  surface  layers  where  all  the 
atmosphere-ocean  interactions  and  energy  transformations  occur.  Since  the 
static  instability  in  the  real  ocean  lasts  only  for  a short  period  of 
time  (hours),  the  simplest  way  of  treating  it  is  an  instantaneous  con- 
vective density  adjustment  mechanism.  Whenever  the  instability  is 
detected  between  layers,  a quasi-homogeneous  mixed  density  of  marginal 
stable  lapse  rate  is  assumed  for  all  unstable  adjacent  layers. 

The  boundary  conditions  for  the  lateral  and  bottom  solid  boundaries 
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are  zero  normal  velocity,  and  zero  normal  heat  or  salt  flux  on  solid 


boundaries  , 


\V  , T ,S 
n n n 


= 0 


(1M 


f 

i i 


where  n denotes  differentiation  along  the  normal  of  the  boundary  surface. 
On  the  equatorial  southern  boundary  in  addition  to  the  no  normal  velocity 
and  no  normal  fluxes  of  heat  and  salt,  the  tangential  free-slip  condition 
is  also  imposed, 


v,u  ,T  ,S  =0  at  4>  = 0. 

’ <t>'  <t>’  <t> 


(15) 


Note  that  the  bottom  topography  and  the  geothermal  heat  flux  have  been 
neglected  at  the  lower  boundary  in  this  preliminary  study.  It  is  known 
that  the  bottom  topography  is  important  to  the  vertically  integrated 
transport  and  can  influence  the  path  of  the  western  boundary  current  in 
the  ocean  (Warren,  1965;  Holland,  1967;  Schulman  and  Niiler,  1970; 
O'Brien,  1971)-  But  we  are  morp  interested  in  the  anomalous  phenomena 
existing  generally  in  the  upper  shallow  layer  of  the  ocean  where  the 
baroclinic  velocity  field  may  not  be  markedly  affected  by  topography. 

At  the  upper  surface  of  the  ocean,  the  "rigid-lid"  approximation 
(Bryan,  1969b)  is  employed  in  order  to  gain  a more  efficient  time- 
dependent  calculation. 

Other  upper  surface  boundary  conditions  are: 


v2(uZ'vz)  = r (t  ,t  } 

o 


(16) 


\ <t> 

where  t , t are  the  longitudinal  and  the  latitudinal  components  of  the 
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surface  stress: 


k2S2  - Sl(E-R) 


(IT) 


where  E - R is  the  net  difference  of  evaporation  and  precipitation;  and 


*2TZ 


PoC 


(18) 


where  Qg  is  the  net  downward  heat  flux  across  the  atmosphere/ocaen  inter- 
face and  C is  the  specific  heat  of  water  at  constant  pressure.  A detailed 
description  of  atmospheric  forcing  function,  namely,  the  wind  stresses  and 
heating  is  in  the  next  section. 

The  boundary  conditions  of  zero  normal  flow  and  zero  heat  and  salinity 
fluxes  through  the  lateral  boundaries,  zero  vertical  velocity  at  the 
balanced  upper  surface  and  at  the  bottom  guarantee  the  satisfaction  of 
integral  constraints  for  mass,  momentum,  and  energy. 

The  wind  stress  and  the  heat  flux.  The  synoptic  net  downward  heat 
flux,  according  to  the  heat  balance  equation  at  the  naviface,  can  be 
written  as 


= Q-r  - Qt.  - Q, 


-*E 


(19) 


Where  Q is  the  net  downward  flux  of  solar  insolation  and  Q , Q , Q are 
I B H T5 

upward  fluxes  of  infrared  radiation,  sensible  heat  and  latent  heat, 
respectively,  as  defined  in  Table  1.  Table  1 also  defines  all  other 
notations  used  in  the  following  formulae. 

The  solar  radiation  flux  can  be  calculated  from  a simplified  formula 
with  empirical  constants  obtained  from  the  atmospheric  climatology 


(Johnson,  et  al. , 1958;  London,  1957;  Vonder-Haar  and  Hanson,  1969)  as 


Qt  = 0.95  Q (0.7^ -0.6N  ) . 
I o c 


The  net  upward  infrared  heat  flux  is  calculated  from 


Q = 0.985  cT  4 (0.59-0. 06  e L/2)  (1-0.6  N 2)  . 

B S a c 


The  sensible  and  latent  heat  are  computed  from 


Qu  = P C C |\V  | (T  -T  ) 
H a H a a s a 


\ * Pa  °E  L IV.I  (V\> 


where  |\V  | is  the  wind  speed  at  10  m above  the  surface.  The  specific 
a 

humidity  is  related  to  the  vapor  pressure  by 


.622 

q = — e 

a 


where  P is  the  mean  atmospheric  surface  pressure.  The  saturated  vapor 
a 

pressure  at  the  ocean  surface  temperature  can  be  calculated  from  the 
Clausius-Clapeyron  equation  (Hess,  1959)* 


L°g10  es  = (9.^051-2553/Tj  . 


The  momentum  flux  at  the  ocjaen-atmosphere  interface  is  expressed  as  the 
surface  wind  stresses, 


(tV)  - pa  cp  |\Val  <Vva) 


where  (u  , v ) are  components  of 'y^.  Following  Deardorff  (1968)  and  using 
data  of  Bussinger,  et  al.  (1971),  the  exchange  coefficients  of  heat,  water, 
and  momentum  are  function  of  atmospheric  stability  as, 


(CD)N  exp( -2  Ri) 


ce  = (ch}n  expt_( 


(CD}N  [1  +^ln  (l'bl  M)] 


CE  (CJpN  ^ + b ln  ^1-b2  Ri^ 
2 


for  stable  cases 


J (Ri  > 


(27a) 


for  unstable  cases 


1 

J (Ri  < 


(27b) 


where  Ri  is  the  bulk  Richardson  number,  a parameter  to  measure  the  sta- 
bility of  the  atmosphere, 

Ri  “ FFi  [<VV  * °-38  Ta  - (28) 

vo  a a 

wherein  p^,  PT,  (CD)N,  b^,  bg,  Z^,  and  are  constants  (Bussinger, 
et  al. , 1971)  as  listed  in  Table  2.  The  coefficient  of  water  particle 
exchange  is  assumed  to  be  the  same  as  that  of  the  exchange  of  heat.  All 
values  of  atmospheric  parameters  are  provided  by  Fleet  Numerical  Weather 
Central  at  Monterey,  California. 

The  hydrostatic  equation.  The  equations  of  state  defines  the  density 
in  terms  of  temperature,  salinity,  and  pressure.  The  departure  of  pres- 
sure at  any  level  in  the  ocean  from  the  vertically  averaged  pressure  is 
uniquely  determined  by  the  density.  Integrating  the  hydrostatic  equation 
downward  with  respect  to  the  vertical  coordinate,  we  obtain 
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p’(z)  = 


Pgdz  - - 


pgdM  dz 


The  vorticity  equation.  The  exclusion  of  kinematic  surface  varia- 
tions makes  the  surface  pressure  impossible  to  calculate.  Instead  we 
compute  the  shear  velocity  components  (u',v')  as  prognostic  variables 
with  the  pressure  terms  substituted  by  the  hydrostatic  equation.  The 
vertical  mean  velocity  (u,v)  is  obtained  from  the  stream  function  by 
solving  the  vorticity  equation.  The  total  velocity  is,  then,  the  combin- 
ation of  the  vertical  mean  velocity  and  the  deviation  from  the  vertical 


mean,  i.e. 


(u, v)  = (u+u’ , v+v' ) . 


The  surface  and  bottom  boundary  conditions  on  w permit  us  to  define 
the  nondivergent  vertically  integrated  stream  function  which  satisfies 
the  continuity  equation. 


A -1 

u = a * 


a -1 

v = a sect. 


where,  for  any  scalar  q. 


(q)  = ~ / (q)dz 

U -D 


is  the  vertically  integrated  mean.  Then  the  predictive  vorticity  equation 
can  be  obtained  by  eliminating  pressure  terms  in  (2),  (3)  as 


sec.  {(^}  „ [(^)cos,]  J . f 


2 
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where  f is  the  Coriolis  parameter  and 


\ uv  ^ 

G = u + — tani)> 
a 


<t>  uv 

G = v - — tan<t> 
a 


The  shear  current.  The  predictive  equations  for  the  vertical  shear 
current  (u’,v'),  from  (2)  and  (3),  subtracting  their  respective  vertical 
means,  become 


, -sec<t>  , „ , \ \ ,a\  a\. 

u!  = p*  + fv'  + F + G - (F  +G  ) 

t p a \ 
o 

. -1  . „ . _<t>  ♦ ,a<1>  a4i. 

v!  = p*  - fu'  + F + G - (F  +G  ) 

t p a <t>  v ' 

o 


where  p’  is  defined  in  (29). 

It  is  obvious  from  (37)  and  (38)  with  the  definition  of  the  vertical 
mean  (33)  that 


(Au;,Av;)  = o. 


With  the  initial  condition  (u',v')  = 0 at  t = 0,  then,  it  is  always  true 


(u',v*)  = 0 


The  prognostic  variables  are  u',  v',  \|r,  T,  and  S predicted  from  (37), 
(58),  (3*0,  (6),  and  (7),  respectively.  The  diagnostic  variables  are  w, 
p,  p obtaining  from  (5),  (8),  and  (29),  respectively. 


5-  Finite  difference  equations 


The  basic  approach  for  arriving  at  the  solutions  is  that  the  system 
of  partial  differential  equations  and  boundary  conditions  in  the  continuum 
are  approximated  by  a set  of  finite-difference  equations  and  bounary  con- 
ditions. As  the  grid  size  approaches  zero,  the  finite  difference  solution 
I 

) obtained  with  any  convergent  scheme  will  approach  the  true  solution. 

However,  since  no  infinitesimal  grid  size  is  used  in  practice,  the  finite 
difference  scheme  chosen  to  compute  the  solution  must  be  capable  of  eco- 
nomically reaching  the  statistical  asymptote  of  the  true  solution.  In 
numerical  simulations  the  governing  equations  are  integrated  beyond  the 
limit  of  mathematical  justification  based  on  local  accuracy.  Therefore, 
it  is  vitally  important  that  certain  integral  constraints  such  as  the 
conservation  of  total  kinetic  energy  during  adiabatic  and  nondissipative 
processes  and  the  conservation  of  total  entropy  and  total  potential  en- 
thalpy during  adiabatic  processes  are  maintained.  In  the  construction 
of  the  model,  we  have  used  the  quadradic  conserving  method  (Arakawa,  1966; 

Bryan,  1969b)  to  maintain  the  conservation  of  integral  constraints. 

I 1 

In  the  NORPAX  model,  the  actual  configuration  of  the  North  Pacific 
Ocean  is  taken  into  account  by  approximating  the  coastline  with  the  mer- 
idians or  the  latitude-circles  passing  through  the  nearest  grid  points 
as  shown  in  Figure  1.  The  origin  of  the  coordinates  is  taken  at  the 
intersection  of  the  mean  sea  surface,  the  equator,  and  the  western  boundary 

of  the  ocean  (X.  = 120°E).  The  whole  North  Pacific  domain  is  from  120°E  to 
o 

80 °W  in  longitude  and  from  0°N  to  65 °N  in  latitude  with  a constant  depth 


■ . , . | . , . • . — „ 


of  4 km.  The  vertical  coordinate  z is  set  for  10  variable-thickness 
layers  as  shown  in  Figure  2.  If  z denotes  the  height,  negative  below 
the  surface,  the  thicknesses  between  the  levels  are  defined  as 


* 


Az 


k+l/2 


= Z. 


k-1 


- z 


V 


2, . . . , K 


(M) 


AZk+l/2  = 2(z  +D) 

IV 


where  D is  the  total  depth.  The  thickness  of  the  kth  level  is  defined 
as 


4\  ' 5 <4\-l/2  + ‘Vl/21’ 


k = 1,...,K 


(^2) 


The  analog  to  (33)  is 
K 

<1  = d Z(q)  AZk  * (U3) 

k=l 

A computational  box  is  one  subvolume  of  the  ocean  with  a finite 
angular  difference  of  A\  in  the  longitudinal  direction  and  A4>  in  the 
meridional  direction  is  shown  in  Figure  3-  We  have  defined 


acosO  , A\ 
«J 


Ay  = aA$ 


(44) 


where  j’  may  refer  to  the  integer  grid  points  or  to  the  half  integer  grid 
points.  Notice  that,  based  on  the  definition  of  (44),  the  x-component  is 
equivalent  to  X-component  and  y-component  is  equivalent  to  ^-component. 


Since  the  computation  is  carried  out  in  flux  form,  the  areas  constituting 


one  subvolume  in  the  ocean  can  be  designated  as 


A1  = AyAzk  ’ 


A2  = Axj'A2k  ’ 


A = Ax  . , Ay 


The  subvolume  under  computation  is 


V,*  ' aXJ'4y4^  • 


For  convenience,  let  us  define 


q 2 1/2  + ^1-1/2^ 


& (q)  “ t(q)x+1y2  - (q) 


where  / may  be  any  of  the  spatial  indices.  For  instance, 


6 ^ ^ 1+1/2,  j+1/2  2 (,?i+l,j+l  + Ti+l,J^  ' 2 (\,j+l  + '?i,j)  ’ 

There  are  five  prognostic  variables,  u’,  v',  f,  T,  S and  three  diag- 
nostic variables,  w,  p,  and  p in  the  model.  In  consistence  with  the  re- 
quirements for  integral  constraints,  the  space-staggered  grids  are  ar- 
ranged in  such  manner  that  T,  S,  and  p are  basically  computed  at  the  in- 
teger grid  points  (i,j,k);  u,v  at  the  half  integer  points  (i+l/2, 
j+l/2,k);  w at  the  (i,j,k+l/2)  point.  The  stream  function  \Jr  is 


two-dimensional  and  is  computed  at  the  (i,j)  points.  All  variables  are 
subscripted  according  to  their  spatial  positions  and  are  superscripted 
with  the  time  step  n,  e.g.,  T11  . The  horizontal  placement  of  vari- 

1,  J 

ables  is  also  shown  in  Figure  3.  All  variables  needed  at  points  other 
than  the  basic  grid  points  during  the  computation  process  are  usually 
defined  by  the  average  values  of  the  parameter  at  the  immediately  avail- 
able points  surrounding  the  temporarily  required  point.  For  example 

_ -j 

Vl/2,j,k  Vl/2,j,k  * 

Momentum  equations.  The  finite  difference  equations  analogous  to 
(37)  for  the  x-component  of  the  vertical  shear  current  with  respect  to 
the  (i+l/2, j+l/2,k)  subvolume  can  be  written  as, 

(u,a)n+1  - 2bAtf(v’a)n+1  = A (^9) 

where  b is  the  implicit  stabilizing  factor  for  the  Coriolis  terms,  t 
the  unit  time  step  and 

A = (u'a)n_1  + 2At(  1-b)  f(v'a)11"1  + 2At  ( Bn+Fn_1+Gn-ft)  (50) 

The  pressure  gradient  force  is  approximated  from  density  alone  according 
to  (29)  as, 

i-.i  \ 

6 - - 6 > r (5D 

o 

where 


kL 


I 


"o,  k = 1 

K-l 


P = 1 V -k 


^k+1/2  g^Zk+l/2  ’ k = 2»***>K 


The  friction  terms,  according  to  (ll)  are, 


Vl  0i,  i 


f»l(6iu)  * i e^A^u)  + |j  (v2Aj6*u) 

! 2tan*  A.Bi(vi)  "1 

. .x  uo  1 


+ v±  ( 1-tan* ) — 

a a cos* 


The  nonlinear  terms,  according  to  (l)  and  (55)  are, 


G = -B^G^u^f]  }-6J{G  [(viA2)J]  )+Bk(uk[wJA  ]")-uv6JA^ 


Note  that  the  right-hand  side  of  (5*0  is,  except  the  last  term,  the  nega- 
tive form  of  the  flux  operator  * The  vertically  averaged  term  is 


0 - J J ( Bn+Fn-1-H3n)Az^  . 


A prognostic  equation  for  the  v' -equation  according  to  (58)  is, 


( v'a)n+1 + 2bAtf(u'a)n+1  = A*  , 


Where  A*  is  similar  to  A,  B*,  F*,  G*  is  similar  to  B,  F,  G,  and  V is 

A 

similar  to  U. 

Thus  the  shear  current  at  all  half  grid  points  ( 1+1/2, j+l/2,k)  are, 


(u’)n+1  = G + 2bAtfQ*  / (l+4b(At)2f2)a 


n+1  p P P 

(V)  = G*  - 2bAtfG  / (144b  (At)  f )a 


21 


Notice  that  in  (57)  and  (58)  the  basic  time  differencing  scheme  used  in 
the  model  is  the  leapfrog  scheme.  This  time  scheme  can  cause  the  tem- 
poral "splitting"  phenomena  in  the  solution  for  a long-term  integration 
process  (Lilly,  19^5) • Therefore,  the  implicit  backward  time  scheme 
(Matsuno,  1966)  is  periodically  used  every  25  time  steps.  The  Coriolis 
term  is  treated  implicitly  which  permits  a time  step  longer  than  the 
inertial  period  (Bryan,  1969b).  For  linear  computational  stability,  a 
stabilizing  factor  b is  chosen  to  provide  a slight  damping  of  inertial 
oscillations,  b = 0.55  is  used  for  the  present  run. 

The  finite  diffence  vorticity  equation  with  respect  to  (i,j)  points 
according  to  (3^)  is 

v2(*D)t  = &1  C~  (6*+£*)]j-&J[^:  (G+^)]i-i  e1!*21*  (59) 


where  G,  F,  G*,  F*  are  the  vertically  averaged  mean  of  g\  f\  g*,  F* 

2± 

equivalently  shown  in  (53)  and  (5M,  6 t indicates  the  derivative  of 
center  difference  with  a separation  of  2Ax,  and 


* ^ 


m 


This  is  a nine-point  finite  difference  Lapacian  operator  which  can  easily 
be  reduced  to  a five-point  Lapacian  operator  if  Ax  = Ay.  With  = 0 

on  all  lateral  boundaries,  (59)  is  solved  by  over-relaxation  in  the  in- 
terior with  an  initial  guess  of  t = 0.  The  in+1  is  advanced  from  \|rn  ^ 
in  pace  with  the  same  time  step  of  the  vertical  shear  current. 

The  vertical  mean  velocity,  with  respect  to  the  computed  points 


(i+l/2,J+l/2)  is. 


A 

U 


A 

V 


= -bJU)Vv 


6l(*)J/Ax 


(61) 

(6l>) 


From  the  continuity  equation  the  vertical  velocity  w can  be  obtained  with 
respect  to  the  (i,j,k)  subvolume  index  as  (i,j,k+l/2), 


Vl/2  * Vl/2  - f [Vi(S;i)  + 


(63) 


where  w.  . = 0 at  the  surface.  Note  that  the  vertical  velocity  is 

i t J > 1/  2 

diagnostically  defined  at  the  top  and  the  bottom  surface  of  the  subvolume 
under  consideration. 

It  is  readily  shown  (Huang,  1973)  that  there  is  no  false  source  or 
sink  of  kinetic  energy  using  the  above  finite  difference  scheme. 

Temperature  and  salinity  equation.  The  temperature  and  salinity  are 
calculated  at  all  integer  grid  points  (i,j,k).  The  finite  difference  form 
of  the  thermal  energy  equation  is 


(Ta)n+1  = (Ta)n_1  + 2At[Qn+(Q'A)n"1] 


where 


Q = -[A181(T1uJ)  + S1(A2Tjv1)] +A5Bk(Tkw) 


(6b) 


(65) 


and 


Q*  = 


k.A 


~ 51(B1T)  + ^ B^A^T]  - A^Ct^ &kT] 


Ax. 


1 2 


5 L&z 


k+l/2 


where  the  eddy  diffusivities  may  be  treated  as  parameters  depending  on  j 
and  k. 
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Similarly,  the  finite  difference  equation  for  salinity  analogous  to 


(7)  has  exactly  the  same  form  if  S is  substituted  for  T.  The  only  source 
and  sink  for  the  potential  energy  in  the  ocean  domain,  besides  the  small 
portion  internally  transformed  from  the  kinetic  energy,  are  the  heat  and 
water  input  through  the  naviface  and  the  diffusive  dissipations. 

Boundary  conditions.  All  the  mesh  grid  points  on  the  solid  lateral 
boundaries,  which  are  approximated  by  the  meridians  and  by  the  latitude 
circles  passing  through  the  nearest  integer  grid  points,  coincide  with 
the  calculating  points  \jr,  T,  and  S.  Understanding  that  the  other  two 
dummy  indices  are  not  specified,  the  zero  flux  condition  at  all  solid 
boundaries  lead  to 


q<+l  ~ qf-l 


where  q is  the  variable  (T, S)  computed  at  the  boundary  and  l is  the  index 
normal  to  the  boundary  at  the  boundary  such  as  i*,  I,  and  J.  However, 


<Jr.  . = 0 


at  all  lateral  boundaries.  The  no-slip  boundary  condition  leads  to 


(u'vW 


However,  for  the  longitudinal  flow  at  the  southern  boundary,  the 
free  slip  boundary  condition  is  also  used.  That  is,  for  j indices, 


(u,v)l/2  - (a,-v)3/2 


2k 


The  downward  fluxes  of  heat,  water,  and  momentum  through  the  naviface  are 


parameterized  as  the  upper  boundary  conditions  indexed  at  k = l/y, 

-5-  6kT  . Hi-  (70) 

“V2  V2  P°° 

~ 5k(S)1/2  - S^(E-R)  , (7D 

5k(u,v)1/2  - (^)  • (7.) 

Hydrostatic  stability.  At  places  where  excessive  upward  heat  flux 
occurs  and  where  evaporation  exceeds  precipitation  by  a large  amount,  a 
layer  of  heavy  cold  and  salty  water  will  be  formed  in  the  upper  layers  of 
the  ocean.  Whenever  the  local  lapse  rate  of  density  is  less  than  the 
vertical  lapse  rate  which  will  ensure  a marginally  stable  vertical  density 
distribution,  the  hydrostatic  instability  exists  in  the  adjacent  layers. 
Since  instability  in  the  real  ocean  usually  lasts  for  a short  period  of 
time,  an  instantaneous  convective  adjustment  mechanism  is  activated  in  the 
model.  Let 


7 - - Q - Pk  ) (73) 

p vVi  /J 

be  the  equivalent  lapse  rate  of  density  between  layers.  A small  positive 

value,  € , is  assumed  since  there  may  exist  a slight  inversion  in  density 

within  the  marginal  stable  limit.  In  the  model  the  lapse  rate  between 

layers  is  computed  and  compared  with  e . If 

P 


•y  > e stable 

P P 


~7  <£ 

P P 


unstable 


where 


P 


is  estimated  from  the  average  density  distribution  in  the  ocean. 


Assuming  the  contribution  to  the  density  variation  from  the  temperature 
and  from  the  salinity  are  equally  important,  the  marginally  stable  cri- 
teria for  temperature  and  salinity  can  be  estimated  from  e and  the  coef- 

P 

ficients  for  temperature  and  salinity  expansions  in  the  mean  state  of 
ocean.  Then,  in  case  instability  is  detected  between  layers,  the  temper- 
ature, salinity,  and  density  in  these  vertical  layers  are  mixed  into 
quasi-homogeneity.  For  instance,  if  instability  exists  in  layer  k-1  and 
k,  the  adjusted  temperature  will  be 


T * = 

k 


T Az,  + (T  - em  Az  / ) Az 
k k v k-1  T k-1/2'  k-1 

Az,  , + Az, 
k-1  k 


T*  = T*  + e Az 
k-1  k e T A\-l/2 


where  the  superscript  * indicates  the  adjusted  value  after  mixing  and 
is  the  marginal  stable  criterion  for  temperature.  These  stability  cri- 
teria and  the  adjusting  processes  are  applied  to  all  layers  in  the  col- 
umn, and  the  procedure  is  repeated  whenever  there  is  more  than  one  un- 
stable subcolumn  until  stability  has  been  reached  for  the  whole  column. 
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4.  Results 


Flexibility  and  efficiency.  In  simulation  studies  using  finite  dif- 
ference schemes  for  long-term  integrations,  not  only  the  conservation  of 
integral  properties  and  the  accuracy  of  solution  in  the  statistic  sense 
are  essential  considerations,  but  the  flexibility  in  choosing  grid  sepa- 
rations and  the  efficiency  in  carrying  out  the  computation  are  also  im- 
portant factors  in  designing  numerical  experiments.  For  the  flexibility 
of  the  model,  we  have  programmed  in  such  a manner  that  the  whole  domain 
is  computed  in  a sequence  of  blocks  which  consist  of  many  computational 
boxes  independent  of  vertical  layers.  We  may  execute  the  experiment  with 
one  version  of  the  model  which  possesses  reasonably  small  space  separa- 
tions, as  small  as  20  km  in  both  longitudinal  and  meridional  directions 
with  as  many  layers  as  needed.  This  is  important  in  the  long  run  for  our 
purpose  of  carrying  out  much-finer-mesh  experiments  in  the  later  state  of 
our  investigation  which  may  involve  the  meso-scale  dynamics.  We  also 
think  our  algorithm  and  computational  sequences  are  optimized  in  such  a 
manner  that  it  took  only  14  sec  of  CDC  7600  CPU  for  one  model  day  with  2.5 
degrees  in  horizontal  and  10  layers  in  vertical.  Depths  of  the  ten  levels 
of  the  2.5-degree  model  are  placed,  respectively,  at  10,  50,  60,  100,  150, 
225,  550,  700,  1500,  and  5000  m below  the  naviface  and  the  constant  maxi- 
mum depth  of  the  ocean  is  assumed  at  4000  m as  shown  in  Figure  2.  It  is 
obvious  that  much  emphasis  has  been  put  on  the  column  of  water  above  the 
permanent  thermocline. 

Three  stages  of  integration.  The  initial  condition  for  the  first  phase 
is  a stable  pattern  of  mean  temperature  and  salinity  with  no  currents.  The 


model  with  5 degrees  in  horizontal  separations  and  7 vertical  layers  was  in- 
tegrated for  more  than  ^0  years  using  8-hour  time  steps  until  a quasi- 
equilibrium state  was  reached.  Figure  4 shows  the  total  kinetic  energy 
fluctuation  at  the  initial  spin-up  stage.  After  10  years  of  integration,  the 
local  rate  of  change  term  was  everywhere  much  smaller  than  other  major  terms 
in  both  the  momentum  and  the  thermohaline  balance  equations,  as  will  be  an- 
alyzed later.  The  results  of  the  model  output  after  1+0  years  of  integration 
was  interpolated  to  2.5  degrees  and  10  layers  to  be  used  as  initial  condi- 
tions for  the  second  phase  of  the  model.  This  version  of  the  model  was  then 
integrated  for  more  than  20  years  until  another  quasi -equilibrium  state  was 
reached.  Then  the  seasonally  varying,  time -dependent  atmospheric  conditions 
were  imposed  as  upper  boundary  forcing  conditions.  The  integration  was  car- 
ried out  for  another  20  years.  The  last  year's  output  was  considered  as  the 
long-term,  seasonally  varying  state  in  the  North  Pacific  Ocean.  In  the  fol- 
lowing papers  (II  and  III),  we  will  describe  in  detail  the  time-dependent, 
seasonal  variability,  and  the  anomaly  experiments. 

Energetic  analysis  under  quasi-equilibrium  state.  The  basic  numer- 
ical scheme  used  in  the  model  was  emphasized  the  conservation  of  integral 
properties  related  to  advection  terms  in  both  the  momentum  and  the  thermo- 
haline equations.  For  each  individual  computational  subvolume,  balances 
of  mass,  momentum,  heat,  and  water  were  properly  taken  into  account. 

After  a long  period  of  time-integration,  a quasi-steady  state  in  the 
ocean  was  reached.  We  will  discuss  a few  representative  points  of  the 
North  Pacific  Ocean  in  depth  to  indicate  the  local  dynamic  balance  of 
all  forces  involved  in  the  system.  Figure  5 shows  the  local  balance  of 


nil  major  physical  processes  included  in  the  formulations,  namely,  the 
Coriolis  term,  the  pressure  gradient  term,  the  horizontal  and  vertical 
viscous  terms,  and  the  horizontal  and  vertical  nonlinear  terms  in  the 
momentum  equation  of  a subvolume  near  the  western  boundary  region. 

Values  indicating  the  relative  magnitudes  of  each  term  in  the  momentum 
equation  have  been  normalized  in  such  a manner  that  the  absolute  sum  of 
all  terms  is  100%.  As  shown  in  Figure  5a  the  dominating  terms  in 
the  momentum  balance  are  the  Coriolis  force,  the  pressure  gradient  force 
and  the  wind  stresses  in  the  surface  layer.  Since  the  x-component  of 
wind  stress  is  much  stronger  than  the  y-component,  the  surface  current 
maintains  an  Ekman-geostrophic  balance  in  the  x-component  of  the  momen- 
tum equations  and  in  near  geostrophic  balance  with  high  frictions  for  the 
y-component  of  the  momentum  equation.  In  the  x-momentum  equation  near 
the  surface  within  the  Kuroshio  region,  the  pressure  gradient  force  and 
wind  stress,  having  relative  weights  about  50%  and  20%,  respectively, 
balance  the  Coriolis  force.  It  is  expected  that  the  nonlinear  effect 
is  relatively  small,  less  than  1%  in  most  places,  as  indicated  in  Figure 
5 because  of  coarse  grids  in  the  model.  The  dominating  flow  in  the 
Kuroshio  region  is  northward  and  there  exists  a strong  shearing  in  cur- 
rents in  the  western  boundary.  As  expected,  the  diffusive  terms  in  the 
y-component  of  momentum  equation  contributes  more  than  10%,  as  shown  in 
Figure  5g,  which  is  much  greater  than  that  in  the  x-component.  Figures 
5b, d show  the  momentum  balance  in  the  lower  layers  of  the  ocean,  at  about 
1000-m  depth.  It  is  obvious  that  the  major  terms  are  in  geostrophic  bal- 
ance. The  summation  of  all  other  terms  contributes  only  less  than  5%. 


, 
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Similar  conclusions  have  been  shown  by  Bryan  and  Cox  (1968)  in  the  vorticity 
balance.  Figure  6 shows  the  momentum  balance  near  the  equatorial  region  at 
2.5 °N,  165 *W.  The  major  forces  are  still  in  Ekman-geostrophic  balance 
in  the  surface  layer  and  in  near  geostrophic  balance  below,  with  the 
horizontal  diffusion  of  momentum  becoming  increasingly  important.  In 
the  x-component  of  momentum  equation,  the  pressure  gradient  force,  which 
is  rather  moderate  in  the  longitudinal  direction,  works  together  with 
the  Coriolis  force  of  balance  x-component  of  the  surface  wind  stress, 
which  is  quite  strong  because  of  the  trade  wind.  The  horizontal  dif- 
fusion of  momentum  as  well  as  the  nonlinear  effect  are  becoming  more 
important  in  the  equatorial  region.  The  former  may  account  for  up  to 
10$  in  the  x-component  and  up  to  20$  in  the  y-component;  and  the  latter 
may  contribute  up  to  5$  at  some  points  ( figures  not  shown) . The  non- 
linear effect  is  small  in  most  of  the  area  in  the  ocean  away  from 
boundaries  but  is  not  negligibly  small  in  the  equatorial  and  in  the 
western  boundary  regions.  The  y-component  of  momentum  equation  has  a 
somewhat  different  balance.  In  the  equatorial  region,  the  y-component 
of  Coriolis  force  is  larger  than  the  x-component  because  the  equatorial 
currents  flow  mostly  in  the  east-west  direction.  The  north-south  den- 
sity gradient  is  larger  than  that  in  the  east-west  direction  while  the 
meridional  component  of  the  trades  is  much  weaker  than  the  zonal  compo- 
nent. The  balance  is  maintained  between  the  Coriolis  force,  working 
together  with  the  meridional  wind  stress  and  the  pressure  gradient  force, 
assisted  by  the  strong  viscous  effect.  In  the  lower  layers  near  the 
equator,  the  x-component  of  the  viscous  force  becomes  important  due 
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mostly  to  large  magnitudes  and  complicated  patterns  of  the  current  system. 
The  viscous  force  together  with  the  Coriolis  force  balance  the  east-west 
pressure  gradient  force.  The  y-componant  of  momentum  in  the  lower  layers 
is  generally  in  geostrophic  balance  as  shown  in  Figure  6d.  In  the  mid- 
latitude open  ocean  about  2000  km  north-west  of  Hawaii  in  the  Pacific 
i wind  drift  current  region,  the  Ekman  dynamics  is  dominating  in  the  x- 

component  while  the  Ekman  geostrophic  balance  is  maintained  in  the  y- 
component  of  momentum  equations  as  shown  in  Figure  7-  At  the  depth  of 
intermediate  water,  there  is  essentially  a geostrophic  balance  with  most 
of  other  processes  being  negligibly  small.  In  summary,  below  the  surface 
layers  geostrophy  is  the  dominating  balance  in  the  present  version  of 
our  model,  while  in  the  surface  layer  the  Ekman  drift  together  with  as- 
sociated convergence  and  divergence  play  an  important  role  in  the  dynamic 
balance. 


The  scale  analysis  shows  that  in  the  density  equation,  consisting  of 
both  the  temperature  and  the  salinity  equations,  the  advective  terms  are 
of  the  same  order  as  other  terms.  It  is  interesting  to  note  that,  though 
the  combined  effect  of  the  horizontal  advection  (V-\VT)  and  the  vertical  ad- 
vection  of  heat  (d/dz(wT))  (or  salinity)  is  comparable  to  the  diffusive 
terms,  the  magnitude  of  the  individual  term  of  the  horizontal  advection 
and  the  vertical  advection  of  heat  (or  salinity)  are  much  lager  than  other 
terms.  Figures  8 and  9 show  the  thermodynamic  balance  in  normalized  rela- 
tive weights  in  percentage  at  two  representative  subvolumes,  one  near  the 
equatorial  region  (near  170°W)  and  the  other  in  the  mid-latitude  open  ocean. 
It  shows  clearly  that  the  most  dominating  process  for  the  local  thermohaline 
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. However,  it  seems 
always  true  that  the  horizontal  advection  terms  and  the  vertical  advec- 
tion  term  are  compensating  each  other  as  they  are  expected  to  be  physi- 
cally meaningful  in  the  continuity  requirement  for  incompessible  fluid 
dynamics.  Divergence  (or  convergence)  will  cause  local  surplus  (or 

1 

deficiency)  of  thermal  energy  which  is  invariably  compensated  by  the  de- 
ficiency (or  surplus)  as  the  result  of  upwelling  (downwelling) . Thus 
the  consequential  net  effect  of  advective  terms  may  still  be  as  small  as 
the  other  processes.  This  general  tendency  holds  all  over  the  ocean 
regardless  of  where  the  computational  subvolume  is.  In  considering  the 
total  balance,  the  surface  heating  turns  out  to  be  as  important  as  the 
advection.  The  former  may  even  be  more  important  than  the  latter  in  the 
open-ocean  region  as  indicated  in  Figure  8 where  the  current  is  rela- 
tively  weak.  From  the  figures  we  notice  that  a real  steady  state  in  the 
density  distribution  has  not  been  reached  though  the  salinity  distribu- 
tion is  almost  steady.  The  term  analysis  of  salinity  ( figures  not  shown) 
shows  the  same  relative  weight  distribution  as  that  of  the  temperature 
distribution. 

The  above  discussion  leads  to  a concerned  question  to  physical 
oceanographers  who  are  seriously  attempting  to  make  in  situ  measurements 
of  heat  content  in  the  ocean.  Since  the  dominating  influence  in  the  local 
thermal  energy  equation  is  the  surface  heating  from  the  atmosphere  and 
the  difference  of  two  large  advective  terms,  small  errors  in  measure- 
ments will  lead  to  a great  uncertainty  in  the  heat  budget  study.  Ac- 
ceptable accuracy  may  be  obtained  with  modern  instruments  and  technology 

I | 


32 


for  the  measurements  of  horizontal  heat  fluxes  and  surface  heating.  But, 
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to  my  knowledge,  the  same  accuracy  for  the  measurement  of  vertical  flux 
of  heat  is  very  difficult  to  obtain  with  the  present  level  of  oceano- 
graphic knowledge. 

Annual  long-term  mean  state.  In  the  initial  phase  of  integration, 

( we  used  a 5-degree  horizontal  separation  with  7 vertical  layers,  of  which 

1 five  were  above  1000  m.  Based  on  the  internal  gravity  wave,  the  Courant- 

Fredericks-Lewy  condition  was  satisfied  in  the  coarse  grid  model  by  tak- 
ing a time  step  of  8 hours.  In  the  second  phase  of  calculation,  the 
model  was  refined  to  2.5  degrees  in  the  horizontal  separations  with  10 
vertical  layers,  of  which  seven  were  above  the  main  thermocline.  A time 
step  of  4.5  hours  was  taken.  Parameters  used  in  the  2.5-degree  version 
of  the  model  are  summarized  in  Table  2.  The  baroclinic  mode  is  predicted 
from  the  prognostic  equations  of  shear  currents  and  the  barotropic  mode 
is  obtained  from  the  stream  function  predicted  by  the  vorticity  equation. 

j 

The  temperature  and  salinity,  hence  the  density,  are  obtained  from  predic- 
tive equations  at  all  time  steps.  The  only  specified  input  to  the  model 
is  the  atmospheric  state  used  in  the  surface  boundary  conditions. 

Mass  transport  stream  function.  Figure  lO  shows  the  contour  plot  of 
the  mass  transport  stream  function  in  a nearly  equilibrium  state  under 
annual  mean  surface  wind  stresses  and  heating  after  more  than  60  years  of 
integration.  Three  large-scale  gyral  circulations  are  well  developed, 
similar  to  the  Munk  solution  (Munk,  1950).  The  anticyclonic  subtropic 
gyre  is  quite  realistically  portrayed  and  the  western  intensification  of 
the  boundary  current  is  clearly  demonstrated.  The  major  port  of  the 


mid-latitude  west  drift  current  constitutes  the  northern  portion  of  the 


subtropic  gyre.  The  California  Current,  the  North  Equatorial  Current  as 
well  as  the  Western  Boundary  Current  are  all  included  in  this  gyral  cur- 
rent system.  The  maximum  total  transport  in  the  Kuroshio  region  near 
Japan  reaches  56  Sverdrups  (56  x 10^  m^/sec)  which  agrees  quantitatively 
with  observed  mean  transport  in  the  area  (Nitani,  1972).  Two  cyclonic 
gyres  are  also  well  simulated,  one  of  which  is  in  the  subartic  region, 
the  other  in  the  tropic  region.  The  subartic  cyclonic  gyre  consists  of 
the  west  drift  current,  the  Alaskan  Current,  The  Aleutian  Current  to- 
gether with  its  extension  through  the  Bering  Sea,  the  Kamchatka  Current, 
and  the  Oyashio  Current.  The  mass  transport  of  the  model  output  in  the 
Kamchatka  Current  region  has  a southward  transport  of  21  Sverdrups  and  in 
the  Oyashio  region,  about  2 Sverdrups.  Reid's  (1973)  analysis  of  the 
Boreas  data  indicates  a geostrophic  transport  across  55 °N  relative  to 
1500  decibars  of  25  Sverdrups  in  winter,  while  calculations  made  from 
summer  data  indicate  only  8 Sverdrups  (Allen,  1964).  In  the  Oyashio  re- 
gion, Hata's  (1965)  data  show  a southward  transport  across  45 °N  of  about 
5 Sverdrups  in  winter  and  about  2 Sverdrups  in  spring  and  summer.  Our 
model  output  agrees  satisfactorily  with  observed  data.  However,  the  sim- 
plification of  configurations  in  the  model  has  ignored  the  existences  of 
islands,  including  the  Aleutian  chain  of  islands,  and  resulted  the  dis- 
appearance of  many  interesting  small-scale  features  in  the  circulation 
pattern.  We  notice  that  the  Aleutian  stream  is  absent  in  the  stream 
function  and  the  Alaskan  gyre  is  flowing  to  the  west  instead  of  to  the 


southwest.  The  mass  transport  in  the  tropic  cyclonic  gyre  with  the  annual 
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mean  atmospheric  forcing  reaches  a maximum  value  of  25  Sverdrups.  The 
gross  feature  of  the  stream  function  patterns  shows  similar  distribution 
to  that  of  the  wind  stress  curl.  Transports  in  the  boundary  regions  be- 
tween the  cyclonic  and  the  anitocyclonic  gyral  circulations  are  weak,  in 
coincidence  with  the  latitude  of  zero  curl  regions  and  the  maximum  trans- 
port occur  at  corresponding  maximum  wind  stress  curl  latitudes  as  pointed 
out  by  Munk  (1950).  The  barotropic  current  deduced  from  the  stream  func- 
tion is  generally  small,  less  than  1 cm/sec  everywhere  in  the  basin  ex- 
cept in  the  western  boundary  region. 

Note  in  Figure  10  there  are  slightly  wavy  streamlines  in  the  mid- 
latitude region  near  Japan.  This  is  probably  because  the  eddy  viscosity 
used  is  still  too  small  for  this  coarse  grid  (Takano,  1976). 

Velocity  patterns.  Figures  11a, b,c  show  current  vector  plots  for  the 

layers  at  a depth  of  10  m,  250  m,  and  2500  m,  respectively.  Figures 

12a, b,c  show  streamline  plots  for  these  same  levels.  In  the  surface  layer 

the  major  swift  currents  occur  in  the  equatorial  region  as  well  as  in  the 

western  boundary  region.  The  maximum  current  is  shown  in  the  eastern 

central  Pacafic  near  the  equator  with  a speed  of  95  cm/sec.  The  Kuroshio 

Current  only  reaches  64  cm/sec  due  mostly  to  the  coarse  grid  meshes.  We 

notice  that  the  western  boundary  current  is  only  barely  resolved  at  the 

horizontal  separation  of  about  175  km.  The  nondimensional  length  scale 

for  a frictional  western  boundary,  L_,  is  proportional  to  the  one-third 

2 1/3 

power  of  the  Ekman  number  (E  = v /2Qa  and  L = E ' ).  The  value  of  L_ 

IF  F 

for  the  second  phase  is  0.012.  The  resolvable  length  scale  for  the  Kuro- 
shio Current  is  340  km  during  the  second  stage  of  computations.  It  is 
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obvious  that  the  grid  separation  is  just  barely  comparable  with  the  obser- 
vational width  of  the  western  boundary  current.  Thus,  the  coarse  grid 
separations  have  led  to  smaller  current  values  than  the  observations, 
while  the  transport  values  and  comparable.  The  scale  width  of  the  iner- 
tial boundary  layer,  h t is  proportional  to  the  square  root  of  the  Rossby 
l/2 

number  (L^.  = Ro  ).  The  nonlinear  effects  were  included  in  the  model 
with  Lj  = 0.0045.  After  comparing  L?  with  L^,  it  is  evident  that  the  in- 
ertial effects  are  smaller  than  the  viscous  effects,  however,  they  are  not 
negligible,  as  ascertained  in  the  previous  momentum  analyses.  The  surface 
currents  in  the  tropics  are  generally  flowing  westward.  The  equatorial 
countercurrent  is  not  clearly  shown,  most  probably  due  to  the  coarse  grid. 
The  Oyashio  merges  into  the  extension  of  the  Kuroshio  to  form  the  North 
Pacific  drifting  current.  The  streamline  plot  for  the  surface  layer  is 
shown  in  Figure  12 c».  It  is  evident  that  the  surface  current  is  dominated 
by  the  wind  pattern.  In  the  equatorial  region,  the  flow  pattern  splits 
into  a northward  branch  to  form  the  Kuroshio  and  a southward  branch  of  the 
equatorial  current  system  near  the  Phillipine  Sea.  The  existence  of  the 
subartic  convergence  zone  is  clearly  shown  along  the  mean  track  of  the 
Kuroshio  and  Oyashio  merging  zone  off  the  western  boundary  at  about  40°N, 
leaning  southward  to  the  eastern  North  Pacific  at  about  35 °N.  A meso- 
scale  anticyclonic  eddy,  recognized  as  the  Eastern  Gyral  by  Sverdrup 
(Sverdrup,  et  al. , 1942)  exists  in  the  mid-latitude  eastern  North  Pacific. 
Figures  lib  and  12b  show  the  current  vector  and  the  streamline  plot  of 
the  layer  at  250-m  depth.  Notice  that  the  current  near  the  equatorial 
region  is  flowing  eastward,  representing  the  combined  flow  of  the  north 
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equatorial  countercurrent  and  the  Cromwell  undercurrent.  The  Kuroshio 

is  well  developed  and  has  the  maximum  velocity  of  55  cm/sec,  the  highest 

value  in  the  layer.  The  250-m  gyre  picture  looks  very  much  like  transport 

maps  of  White  (1975).  Note  that  the  surface  convergence  line  exists  at 

the  center  of  the  geostrophic  gyre  at  250  m and  that  the  eastern  surface 
l 

i gyre  does  not  show  up  at  depth.  The  Kuroshio  Current,  the  west  wind 

drift,  the  California  Current,  and  the  North  Equatorial  Current  are  pre- 
dominantly and  easily  identified  with  the  anticyclonic  gyre.  The  tropic 
cyclonic  gyre  consisting  of  the  North  Equatorial  Current,  the  equa- 
torial countercurrent  together  with  the  equatorial  undercurrent,  is 
clearly  demonstrated.  Other  less  obvious  currents,  such  as  the  Oyashio 
and  the  Alaska  Currents,  are  indicated  in  the  circulation  pattern.  Figures 
11a, b and  12a, b simulate  the  gross  nature  of  the  circulation  patterns  in 
the  upper  layers  of  the  North  Pacific  Ocean  as,  for  example,  shown  in 
Figure  15  which  indicates  the  surface  mean  current  in  more  detail.  The 
simulated  current  near  700  m is  rather  weak  and  the  circulation  indi- 
cates a different  pattern  ( figure  not  shown) . In  the  deep  layers  cur- 
rents are  slightly  stronger  than  in  layer  7 and  the  circulation  pattern 
below  the  depth  of  1500  m,  as  shown  in  Figures  11c  and  12c,  is  opposite 
to  that  of  the  upper  layers.  The  current  beneath  the  Kuroshio  is  flowing 
southward  at  a much  slower  speed  of  4 cm/sec  and  there  is  clearly  a weak 
cyclonic  circulation  in  the  subtropic  region.  Deep  water  is  slowly  flow- 
ing southward  from  the  high  latitudes  in  the  western  North  Pacific  and 
northward  in  the  eastern  North  Pacific,  and  that  generally  agrees  with 
the  deep  circulation  pattern  deduced  from  the  geostrophic  calculation 
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of  Reid  (1976,  personal  communication). 

Figure  l4  shows  the  vertical  velocity  contours  at  the  bottoms  of 
aforementioned  three  layers  (i.e.,  at  20  m,  300  m,  and  2200  m depths, 
respectively)  in  the  model  ocean.  The  solid  contour  indicates  upwelling 
and  the  dotted  line,  downwelling.  There  are,  in  general,  strong  upwell- 
ings  along  the  equator,  especially  in  the  eastern  part  of  the  mid-tropic 
region  where  the  maximum  upwelling  reaches  5 m/day.  Upwellings  in  the 
mid-latitude,  western  boundary  region  are  relatively  weak  in  the  surface 
layer.  There  also  exist  weak  upwellings  along  the  California  coast.  Weak 
downdwellings  exist  in  most  subtropic  and  eastern  tropic  regions  away  from 
the  equator.  This  pattern  of  equatorial  upwellings  and  subtropic  down- 
dwellings  extends  to  a depth  of  more  than  700  m.  In  deeper  layers,  up- 
wellings along  the  western  boundary  and  downwellings  in  the  high  latitudes 
and  in  the  north-western  part  of  the  North  Pacific  become  even  stronger 
at  the  depth  of  1500  m.  The  maximum  upwelling  reached  3*5  m/day  in  the 
western  boundary  at  a depth  of  200  m and  the  maximum  downwelling  of  ap- 
proximately 1 m/day  at  the  high  latitude  region  occurred  at  a depth  of 
530  m.  The  equatorial  upwelling  decreases  below  200  m and  a slight  up- 
welling was  indicated  in  the  tropics  at  the  1000-m  depth.  However,  below 
1000  m,  the  general  trend  of  upwelling  and  downwellings  generally  changed 
signs.  There  exist  moderate  downwellings  in  the  western  boundary  at 
1000  m and  weak  upwellings  in  the  high  latitudes.  At  the  depth  of  2500  m, 
weak  upwellings  become  evident  in  most  of  the  northern  boundaries  and 
moderate  downwellings  exist  in  the  equatorial  region.  There  are  weak 
but  broad  upwellings  in  the  interior  of  the  tropics.  From  the  previous 
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figures  we  notice  that  major  large-scale  features  in  the  model  match 
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nicely  with  observed  mean  current  patterns  and  upwelling  distributions 
in  the  North  Pacific  Ocean. 

Density  field.  Density  is  a diagnostic  parameter  obtained  from  the 
temperature  and  salinity  according  to  the  equation  of  state.  Several 
formulations  based  on  empirical  data  are  available  (e.g.,  Eckart,  1958; 
Bryan  and  Cox,  1970;  Friederick  and  Levitus,  1970).  In  the  present  run 
the  simplest  form  was  used, 

p = p0[l-a(T-To)  + r(S-So) ] (76) 

where  values  of  To,  So,  a,  and  7 are  chosen  from  appropriate  data  as 
shown  in  Table  2.  Since  both  T and  S are  prognostic  variables  in 
the  model,  we  have  examined  the  temperature  and  salinity  field  instead  of 
the  density  field.  Figure  15  shows  the  temperature  patterns  at  depths 
of  10  ra,  250  m,  and  5000  m,  respectively.  The  surface  temperature  dis- 
tribution (Figure  15a)  shows  an  equatorial  cool  region  extending  from  the 
east  boundary  to  l80°,  corresponding  to  the  equatorial  upwelling  region. 
The  broadest  cold  water  extends  to  5®N  with  a minimum  of  23 °C  at  the 
equator.  There  is  a broad  warm  water  region  in  the  western  tropic  and  a 
narrow  one  near  the  eastern  boundary  north  of  the  equator.  A sharp  merid- 
ional temperature  gradient  exists  in  the  mid-latitude  off  the  western 
boundary  along  55 °N.  The  upwelling  along  the  equator  does  not  penetrate 
deeply  around  500  m,  warm  water  masses  in  the  western  and  eastern  tropics 
can  still  clearly  be  discerned  and  a sharp  temperature  gradient  along  the 
western  boundary  develops.  At  1000  m (figure  not  shown),  the  temperature 
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is  quasi -homogeneous  and  is  slightly  warm  in  the  south.  The  deep  water 


at  3000  m shown  in  Figure  15c  indicates  a spatial  homogeneity  with  a 
slight  cooling  in  the  north.  Because  of  our  limited  knowledge  concern- 
ing the  interior  of  the  ocean,  verifications  for  the  model  ouput  are 
rather  difficult.  Nevertheless,  certain  climatological  atlases  are 
available  for  comparisons.  In  general,  the  distribution  agrees  well  with 
the  mean  sea  surface  temperature  pattern  in  the  North  Pacific  Ocean  as 
shown  in  Figure  1 6 (Masuzawa,  1972).  In  the  mid-latitude  and  tropic  re- 
gionSj  the  agreement  between  Figures  15a  and  l6  is  very  good.  However,  a 
close  examination  reveals  that  the  observed  winter  temperature  at  high 
latitudes  is  slightly  lower  than  the  model  temperature,  which  is  much 
closer  to  the  summer  values  (Figure  l6a) . Since  the  temperature  changes 
depend  mostly  on  the  constrained  atmospheric  conditions,  especially  the 
air  temperature,  the  annual  mean  atmospheric  data  used  in  the  model  as 
upper  boundary  conditions  is  probably  biased  to  the  summer  conditions 
mostly  because  ther  were  more  observations  during  the  summer  season.  The 
similar  temperature  patterns  extend  down  below  to  100  m with  a slight 
overall  decrease  of  about  2°C  at  the  60-m  depth  and  about  7°C  at  the  100-m 
depth,  except  the  cool  water  along  the  equator,  that  has  been  narrowed 
to  1 or  2 degrees  in  latitude  at  the  60-m  depth  and  disappeared  at  100  m. 

Temperature  distributions  for  longitudinal  vertical  cross  sections 
at  5°N,  15 °N,  25 °N,  and  35 °N  are  shown  in  Figure  17.  There  is  always  a 
sharp  vertical  temperature  gradient  about  0.5°C/m  above  the  main  thermo-  < 
cline,  which  is  about  300  m in  the  tropics.  It  is  about  0.015°C/m  in  the 
intermediate  water  from  the  main  thermocline  down  to  the  depth  of  1000  m. 
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The  temperature  is  rather  homogeneous  with  a weak  vertical  gradient  about 
0.0007 °C/m  below  the  depth  of  1000  m.  Water  is  generally  warm  in  the 
west  North  Pacific  (32°C  to  3^°C),  cool  in  the  middle-east  North  Pacific 
(about  2h *C  to  25 °C),  and  becomes  warm  again  near  the  eastern  boundary 
(about  29*C  to  30°C)  in  the  surface  layer  at  5°N  near  the  topics.  The 
eastern  warm  mass  extends  westward  to  135 °N  and  more.  The  surface 
mixed  layer  is  about  50  m to  100  m deep  in  the  mid-latitude  as  shown  in 
Figure  17d.  The  observed  corresponding  temperature  distributions  in  the 
longitudinal  cross  sections  of  Figure  17  are  shown  in  Figure  18,  obtained 
from  Moromtsev  (1963)  and  from  Masuzawa  (1972).  The  general  feature  of 
the  model  output  agrees  qualitatively  well  with  observed  and  atlas  data. 
The  same  order  of  magnitude  of  vertical  temperature  gradients  exists  in 
the  upper  350  m.  However,  there  are  no  spatially  wavy  distributions  of 
temperature  and  no  double  thermoclines  in  the  simulated  temperature  field 
as  have  been  shown  in  the  atlas.  This  is  probably  due  to  the  coarse 
vertical  separations  in  the  model.  The  model  of  the  eastern  North  Pa- 
cific produced  more  warm  water  in  the  middle  latitude,  where  the  iso- 
therms rise  very  gently  to  the  east  in  the  intermediate  water  depth 
(300  m to  700  m),  as  indicated  in  Figures  17c, d,  than  that  of  the  atlas 
where  the  isotherms  rise  rather  sharply  toward  the  east  as  shown  in 
Figures  l8c,d. 

Four  meridional  cross  sections  of  temperature  distribution  of  the 
model  output  are  shewn  in  Figure  19  and  four  corresponding  cross  sections 
of  oceanographic  surveyed  profiles  are  shown  in  Figure  20.  The  gross  fea- 
tures of  the  two  figures  bear  general  agreement,  except  the  isotherms  in 


the  model  output  are  much  smoother  than  the  observed  ones.  The  major 


discrepancy  occurs  in  the  high  latitude  region  above  1)5°N  where  the 
simulated  temperature  is  warmer.  The  1|  °C  isotherm  never  reaches  the 
surface  in  the  model  as  it  does  in  the  observed  data.  This  is  mostly 
caused  by  the  annual  mean  atmospheric  forcing  conditions  which  are 
biased  by  northern  summer  observations.  Warm  air  temperatures  in  high 
latitudes  used  as  constrained  boundary  conditions  result  in  a corre- 
spondingly warm  oceanic  temperature  distribution  in  the  subartic  re- 
gion. The  overadjustment  of  the  convective  mechanism  in  the  unstable 
water  column  also  affects  the  mixed  layer  depth  which  is  rather  deep  in 
the  model. 

Simulated  salinity  distributions  for  three  layers  at  surface,  350  m, 
and  near  the  bottom  (3000  m)  are  shown  in  Figures  21a, b,c,  respectively. 

The  surface  salinity  shows  an  elongated  maximum  zone  above  35°/°°  between 
20 °N  and  30 °N  and  declines  polarward  and  equatorward  as  indicated  in  Figure 
21a.  In  the  eastern  tropic,  there  is  also  a high  salinity  tongue  near  the 
boundary  and  extending  westward.  The  north-south  homogeneity  is  indicated 
in  the  western  boundary  region.  At  the  depth  of  intermediate  water  level, 
salinity  values  are  slightly  higher,  around  31+-l+°/oo  to  35-2°/oo,  than  that 
of  the  surface  layer  but  with  similar  patterns  as  shown  in  Figure  21b.  The 
salinity  in  the  bottom  water  is  rather  homogeneous  with  a slight  high  in  the 
northwest.  Reid's  (1969)  surface  salinity  distributions  are  shown  in  Figure  22. 
General  agreements  in  large-scale  features  between  Figures  21a  and  22  are  satisfac- 
tory. Three  meridional  vertical  cross  sections  at  l80°,  105 °W,  and  l65°E  are  shown  in 


Figures  23a, b,c,  respectively,  for  the  simulation  and  in  Figures  2Ua,b,c, 


for  the  observed  atlas  data  (Moromtzev,  1963)*  The  gross  nature  in  the 
corresponding  cross  sections  in  Figure  23  and  Figure  24  agrees  qualita- 
tively well.  In  the  western  North  Pacific  along  105 °Wf  a maximum  near 
20°N  and  a minimum  near  5°N  are  shown  in  Figures  23a  and  24a.  Along  the 
l80°  cross  section,  a large  amount  of  low-salinity  but  high-density  water 
sinking  near  40°N  to  an  intermediate  water  depth  around  500  m to  1000  m 
is  clearly  indicated  in  both  Figures  23b  and  24b.  There  is  a shallow 
salinity  minimum  in  the  vertical  at  the  depth  around  100  m,  a maximum  at 
the  depth  around  300  m and  a prominent  minimum  at  the  intermediate  water 
depth.  These  features  are  also  clearly  shown  in  the  eastern  central 
North  Pacific  in  both  Figures  23c  and  24c.  Certain  discrepancies  do 
exist  between  the  simulated  and  the  observed  salinity  distributions. 

The  simulated  salinity  in  the  high-latitude  region  north  of  40°N  is 
about  l°/o°  higher  than  atlas  data  and  the  maximum  salinity  zone  is  closer 
to  the  eastern  boundary  than  is  the  observed  one.  This  is  mostly  caused 
by  the  inaccuracy  of  the  evaporation  and  precipitation  data  imposed  as 
the  upper  boundary  condition.  These  data  were  obained  from  Jacobs  (1951). 
Recent  studies  (e.g.,  Reed  and  Elliot,  1973)  have  indicated  that  much 
high  precipitation  occurs  in  the  subartic  region. 

Pressure  field  and  dynamic  topography.  The  surface  pressure  can  be 
calculated  explicitly  from  the  vertically  integrated  momentum  equations, 
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where  PQ  is  the  surface  pressure,  ^ is  the  stream  function  defined  in 

(31)  and  (32).  F and  G are  defined  in  (9),  (10),  (31*,  (55).  The  p can 

o 

be  obtained  by  integrations  with  an  arbitrary  constant  appropriately  ad- 
justed with  a known  pressure  point.  Then  the  pressure  at  all  depths  can 
be  found  by  integrating  the  hydrostatic  Eq.  (29).  The  surface  barotropic 
pressure  above  the  balanced  level  surface  can  be  approximated,  excluding 
the  very  small  surface  baroclinic  and  kinematic  effects,  as 


Po  = Pa  + P0gTl  (79) 

where  q is  the  surface  elevation  and  Pa  is  the  atmospheric  surface  pres- 
sure. 

The  pressure  field  at  the  surface  (10  m)  and  at  the  depths  of  225  m 
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and  3000  m are  shown  in  Figure  25  in  units  of  0.1  m /sec  . The  pres- 
sure field,  Figure  25a,  shows  sharp  horizontal  gradients  in  the  Kuroshio 
region,  indicating  a strong  northward  current.  The  pressure  gradients 
outside  of  the  boundary  indicate  a broad  eastward  flow  in  the  mid-latitudes 
and  a weak  south-westward  flow  in  the  tropics.  At  the  depth  of  225  m,  the 
basic  pattern  of  surface  flow  still  exists  at  this  depth  except  with  more 
northward  flow  in  the  subartic  and  more  westward  flow  in  the  tropics.  Near 
the  bottom  layer  at  a depth  of  3000  m,  the  geostrophic  circulation  is 
opposite  of  that  in  the  upper  layers.  There  is  southward  flow  beneath 
the  Kuroshio  Current  and  a broad  westward  flow  under  the  North  Equator 
Current. 


The  dynamic  topography  of  the  sea  surface  based  on  observational 
density  field  with  respect  to  a 1000-db  is  shown  by  Wyrtki  (197^)  and  the 


mean  dynamic  height  in  the  North  Pacific  Ocean  is  reproduced  as  shown  in 
Figure  26.  Though  detailed  information  is  not  portrayed  in  the  computed 
pressure  field,  especially  near  the  equator,  general  features  of  the  sur- 
face pressure.  Figure  25a  as  the  relative  field  of  the  surface  elevation 
compares  very  favorably  with  annual  mean  dynamic  topography  in  the  North 
Pacific  as  shown  in  Figure  26.  Note  that  the  dynamic  height  in  Figure  26 
corresponds  to  the  surface  pressure  values  in  Figure  25  divided  by  the 
gravitational  acceleration,  and  the  difference  is  just  an  arbitrary  inte- 
gration constant. 

Energetics . The  energy  conservation  for  the  model  has  been  dis- 
cussed briefly  in  the  last  paragraph  of  the  formulation.  The  rate  of 
change  of  kinetic  energy  per  unit  volume  can  be  expressed  symbolically 
as 

K = K ' + K = <G>  + <p~K>  + <t>  + <DV>  +•  <DW>  (80) 

where  K and  K’  are  the  rate  of  change  of  vertical  mean  and  deviation  from 
the  mean  kinetic  energy  per  unit  volume;  G denotes  the  nonlinear  effect; 
(pr~K),  the  transformation  between  the  potential  energy  and  the  kinetic 
energy  due  to  buoyant  effect;  t,  the  work  done  on  the  surface  of  the  ocean 
by  wind  stress;  (DV),  the  dissipation  of  kinetic  energy  due  to  horizontal 
eddy  viscosity;  and  (DW),  the  dissipation  of  kinetic  energy  due  to  vertical 
eddy  viscosity.  The  rate  of  change  of  potential  energy  per  unit  volume  can 
be  expressed  as 


<GE>  - <p~K>  + <DP> 


(81) 


where  P = pg(z+D)  is  the  potential  energy  defined  as  before;  (GP)  denotes 
the  advection  of  potential  energy;  (pr-K)  the  conversion  of  potential 
energy  to  kinetic  energy  or  vice  versa;  and  (DP),  the  diffusion  of  poten- 
tial energy  including  all  thermocline  sources  and  sinks.  Note  that  if 
integration  is  carried  out  for  the  whole  mechanically  closed  domain,  the 
advection  effects  such  as  G and  (GP)  contribute  no  work.  However,  redis- 
tribution of  kinetic  energy  between  the  baroclinic  mode  K'  and  the  baro- 
tropic  mode  K does  happen  because  of  the  nonlinearity.  Note  also  that 
the  conversion  between  the  potential  and  the  kinetic  energy  occurs  only 
in  the  baroclinic  mode  because  the  ocean  bottom  is  flat  in  the  present 
study  (Wiin-Nielsen,  1962;  Holland,  1973). 

The  energetic  diagram  showing  the  magnitude  and  direction  of  work 
done  per  unit  volume  in  ergs  cm  ^ day  , together  with  the  total  poten- 
tial and  the  total  barotropic  and  baroclinic  energy  in  erg/crr?  are  shown 
in  Figure  27.  Potential  energy  is  computed  with  reference  to  the  ocean 
bottom.  It  is  clear  from  Figure  27  that  the  kinetic  energy,  from  both 
the  barotropic  and  the  baroclinic  modes,  has  reached  a steady  state  but 
the  potential  energy  is  still  losing  energy  at  a rate  of  2.5 
ergs/cm^/day.  In  the  North  Pacific  model,  a major  portion  of  the  surface 
area  is  in  the  tropic  and  subtropic  region.  This  leads  to  a relatively 
large  amount  of  input  heat  and  output  evaporation  across  the  ocean  sur- 
face as  a whole.  These  net  thermal  and  haline  effects  dominate  the 
equilibrium  potential  energy  state,  the  former  decreases  the  potential 
energy  while  the  latter  increases  it.  As  shown  in  the  energy  diagram, 
the  thermal  effect  is  very  influential.  A small  portion  of  potential 
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energy  is  increased  due  to  the  net  upwelling  induced  by  the  mechanical 
input  of  wind  stresses.  As  pointed  out  by  Haney  (197^),  consistent  with 
Holland's  (1972)  result,  the  energy  transformation  is  from  the  kinetic  to 
the  potential  energy  when  the  thermal  boundary  layer  is  greater  than  the 
viscous  boundary  layer  as  it  is  in  the  case  of  the  present  model,  based 
on  scale  considerations.  The  diffusive  effects  also  act  on  the  positive 
j side  of  the  potential  energy.  Since  the  model  has  reached  a quasi - 

equilibrium  state,  the  energy  loss  due  to  the  convective  adjustments  under 
unstable  conditions  is  small.  The  kinetic  energy  is  derived  from  the 
surface  wind.  Most  of  the  mechanical  energy  input  goes  to  the  baroclinic 
mode.  There  is  a small  amount  of  energy  conversion  between  the  baro- 
clinic mode  and  the  barotropic  mode  due  to  the  nonlinear  effect. 

Energy  transport.  One  interesting  and  important  problem  in  the 
numerical  ocean  simulation  is  to  find  out  how  much  energy  is  transported 
by  the  ocean  in  the  atmosphere/ocean  system.  In  general,  thermal  energy 
must  be  transported  poleward  from  low  latitudes  to  high  latitudes.  The 
total  northward  transport  of  heat  across  a latitudinal  circle  in  the 
ocean  domain  can  be  expressed  as 
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where  terms  in  the  square  bracket  are  the  barotropic  and  baroclinic  part 


of  the  meridional  velocity  and  the  north-south  diffusion  due  to  tempera- 
ture gradient  and  and  X^.  are  the  longitudes  of  the  western  and  the 
eastern  boundaries,  respectively.  The  results  are  shown  in  Figure  28. 


is  transported  by  the  baroclinic  advection  (Hv')  which  has  approximately 


the  same  magnitude  as  the  energy  flux  transported  by  the  total  current 
( Hv) . The  diffusion  of  heat  (Hd)  contributes  relatively  little  in  the  low 
and  high  latitudes  and  is  equatorward  in  the  tropics.  However,  in  the 
mid-latitude  region  where  the  meridional  temperature  gradient  is  large 
and  the  current  is  mostly  zonal,  the  transport  of  heat  energy  by  dif- 
fusive processes  is  important.  The  total  rate  of  the  northward  thermal 
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energy  transport  is  about  2 x 10  cal/sec.  Maximum  transports  occur  at 
latitudes  between  20  ®N  and  35 °N. 

Since  the  poleward  energy  transport  by  ocean  currents  is  important 
in  the  dynamics  of  climate  fluctuations,  it  has  been  studied  for  many 
years  by  many  scientists  (e.g.,  Jung,  1952;  Sverdrup,  1957;  Budyko,  1958; 
Albrecht,  I960;  Bryan,  1962).  Vonder-Haar  and  Oort  (1975)  have  made  a 
new  estimate  of  the  annual  poleward  energy  transport  based  on  recent 
measurements  of  the  earth's  radiation  budget  from  satellites,  together 
with  atmospheric  energy  transport  summaries  based  on  rawinsonde  data. 


Their  results  are  shown  in  Figure  29*  The  peak  oceanic  heat  transport 
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(2.5  x 10  cal/year)  occurs  near  20°N  and  accounts  for  7 4$  of  the  energy 
transport  at  that  latitude.  On  the  average,  the  ocean  contributes  about 
bQrfo  in  the  northern  hemisphere  for  the  required  energy  in  the  atmosphere/ 
ocean/earth  system.  Sellers'  (1965)  estimation  based  on  the  surface  energy 
budget  of  Budyko  (1963)  is  also  shown  in  Figure  29,  as  BS.  As  pointed  out 
by  Vonder-Haar  and  Oort,  Sellers'  estimation  is  relatively  lower.  The  com- 
puted energy  transport  in  the  North  Pacific  Ocean  model,  indicated  as 
NPO,  is  also  shown  in  Figure  29.  Though  quantitative  comparison  between 


the  observational  estimates  and  the  model  output  without  the  energy  trans- 
port from  the  Atlantic  and  the  Indian  Ocean  is  not  meaningful,  it  is  certain 
that  the  total  simulated  oceanic  energy  transport  including  the  above- 
mentioned  transports  would  be  closer  to  the  recent  estimate  by  Vonder- 
Haar  and  Oort  (1973). 

! 
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5.  Summary  and  conclusions 


Based  on  Haney's  (197*0  idealized  basin  model,  a North  Pacific  Ocean 
model  has  been  developed  for  the  study  of  the  physical  nature  and  dynamic 
behaviors  of  large-scale,  long-term  fluctuations  in  the  North  Pacific 
Ocean.  During  the  first  stage  of  time  integration,  a coarse  grid  mesh 
with  horizontal  separation  of  5 degrees  and  7 layers,  was  used  to  obtain 
the  long-term  mean  state  of  the  circulation.  In  the  second  stage,  the 
model  was  refined  to  2.5  degrees  in  horizontal  separations  with  10  layers 
in  vertical.  Annual  mean  observed  atmospheric  data  were  imposed  as  the 
upper  boundary  conditions.  After  about  60  years  of  integration  in  the 
second  phase,  a quasi-equilibrium  state  in  the  ocean  was  reached  and  the 
time -dependent,  seasonal  atmospheric  conditions  were  imposed  to  obtain 
seasonal  responses  and  time -dependent  variations  in  the  ocean.  Another 
20  years  or  more  integration  was  carried  out  and  daily  values  of  the  last 
years  were  used  as  the  mean  seasonal  state.  Prediction  experiments  were 
carried  out  using  the  observed  anomalies  together  with  the  mean  seasonal 
state  in  the  ocean  for  the  understanding  of  the  dynamics  involved  in  the 
generation,  propagation,  evolution,  and  dissipation  of  thermal  anomalies 
which  were  observed  in  the  North  Pacific  Ocean.  The  present  paper  deals 
with  the  long-term  mean  state  in  the  completion  of  the  second  phase. 

Dynamic  analyses  show  that  geostrophic  balance  is  approximately  main- 
tained all  over  the  ocean  in  the  lower  layers.  In  the  surface  layer,  the 
momentum  balances  are  in  general  between  the  Coriolis  force  and  the  sur- 
face stress  in  the  zonal  equation  and  between  the  Coriolis  force,  the 
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surface  stress,  and  the  pressure  gradient  force  in  the  meridional  equa- 


tion. Internal  vertical  eddy  viscosity  contributes  about  10$,  to  the 
balance.  The  nonlinear  terms  are  generally  small  but  not  negligible. 

In  the  heat  and  salinity  equations,  the  convergency  of  the  advective 
fluxes  of  heat  and  salinity  are  the  dominating  terms.  This  leads  to  the 
question  as  to  how  accurately  the  energy  budget  in  a measurable  volume  of 
water  in  the  real  ocean  can  be  calculated.  Since  ohe  horizontal  and  ver- 
tical fluxes  tend  to  balance,  small  errors  involved  in  the  measurement  of 
the  transports,  especially  the  vertical  one,  may  result  in  large  residual 
errors  in  the  heat  budget. 

The  vertically  integrated  stream  function  shows  three  gyre  circula- 
tors in  the  North  Pacific  Ocean,  one  anticyclonic  gyre  in  the  subtropic, 
and  two  cyclonic  gyres  in  the  tropic  and  subartic  regions,  respectively, 
similar  to  the  solutions  of  a diffusive  ocean  (Munk,  1950).  Layer-by- 
layer  velocity  patterns  show  surface  flow  toward  the  west  and  deep  flow 
toward  the  east  near  the  equator,  indicating  the  simulation  of  the  South 
Equatorial  Current  and  the  Cromwell  Current.  There  is  a northward  Kuro- 
shio  Current  above  700  m and  a returning  southward  flow  underneath  it. 

The  layer-by-layer  streamline  patterns  show  large-scale  anticyclonic  cir- 
culations above  1000  m and  opposite  patterns  below.  There  are  strong  up- 
wellings  near  the  equator  and  at  the  western  boundary  in  the  upper  two 
layers  and  downwellings  in  the  lower  layers  in  the  equatorial  region.  The 
simulated  gross  nature  of  temperature  and  salinity  distributions  compare 
well  with  observed  field  atlas  data.  The  surface  pressure  field  also  com- 
pares satisfactorily  with  dynamic  height  distribution  based  on  density 
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observations . The  northward  transport  of  thermal  energy  is  accomplished 


mostly  by  the  baroclinic  current  and  the  total  thermal  energy  transport  is 
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rather  large  (1*1x10  cal/year).  However,  both  temperature  and  salinity 
values  from  the  model  are  higher  than  observed  mean  data  in  high  latitude. 
This  is  mostly  due  to  the  higher  annual  mean  air  temperature  and  lower 
evaporation  and  precipitation  values  imposed  upon  the  model  than  these  in 
) the  real  atmospheric  field. 

The  present  regional  North  Pacific  model  demonstrated  its  ability  to 
portray  the  large-scale  mean  circulation  patterns  and  mean  density  field 
in  the  North  Pacific.  The  major  objective  of  performing  the  anomaly  pre- 
diction experiments  for  the  understanding  of  anomaly  dynamics  in  the  ocean 
justifies  the  step-by-step  approach  to  obtain  the  time-dependent,  long- 
term mean  ocanic  state.  A series  of  numerical  studies,  imposing  anom- 
alous heat  fluxes  and  varying  wind  stresses  as  well  as  the  anomalous 
initial  state  in  the  ocean  have  been  carried  out.  The  follow-up  papers, 

II  and  III,  will  deal  with  time -dependent  oceanic  state's  and  anomaly 
experiments. 
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Table  1 


Description  of  symbols  in  ocean  model 


time 

longitude 

latitude 

height 

angular  speed  of  the  earth 

radius  of  the  earth 

acceleration  of  gravity 

zonal  component  of  the  current 

meridional  component  of  the  current 

reference  horizontal  velocity 

vertical  component  of  the  current 

temperature 

salinity 

density 

eastward  frictional  force  per  unit  volume 
northward  frictional  force  per  unit  volume 
horizontal  gradient  operator 

step  function  incorporated  with  density  stability 
nonlinear  component  in  the  zonal  direction 
nonlinear  component  in  the  meridional  direction 
downward  heat  flux 
incoming  solar  radiation 
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Table  1 (concluded) 


B 
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e 
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c 
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V 

a 
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p 
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C ,C  ,C 
D’  E H 


Ri 

X 


JP 


back  radiation 
sensible  heat 
latent  heat 

incoming  solar  radiation  with  a clear  sky- 

cloudiness  in  tenths 

vapor  pressure  of  air 

saturated  vapor  pressure 

ocean  surface  temperature 

air  temperature 

Stefan- Boltzmann  constant 

latent  heat  of  evaporation 

air  density 

wind  velocity 

specific  humidity 

air  pressure 

drag  coefficient  for  momentum,  latent  heat,  and  sensible  heat, 

respectively 

bulk  Richardson  number 

longitudinal  component  of  wind  stress 

latitudinal  component  of  wind  stress 
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Values  of  parameters  and  constants  in  the  ocean  model 


Table  2 (continued) 


constants  in  drag  coefficients  b 53-2  ergs 


Table  2 (concluded) 


List  of  figures 
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‘ 
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Figure  1. 

Figure  2. 

Figure  3* 
Figure  it. 
Figure  5- 


Figure  6. 
Figure  7- 

Figure  8. 


Figure  9« 
Figure  10. 

Figure  11. 


NORPAX  model  2.3-degree  x 2.3-degree  grid  mesh  and  param- 
eters computing  points:  x temperature  and  salinity  points,  • 
velocity  points. 

Level  depths  for  vertical  layers  and  layer  thicknesses  at  the 
integer  k and  the  half  integer  k+l/2  level. 

A typical  computational  box,  i.e.,  a subvolume. 

Kinetic  energy  vs.  time  plot  in  the  spin-up  stage. 

Relative  weight  of  terms  (%)  for  the  momentum  equation  in 
the  western  boundary  regions:  a and  b the  x-component  of 
momentum  balance  for  the  layers  near  10 -m  and  1000-m  depths; 
c and  d the  y-component  of  momentum  at  the  same  depths. 

Same  as  Figure  5 except  for  the  equatorial  region  near  170°W. 
Same  as  Figure  5 except  for  the  open  ocean  in  the  middle 
latitude. 

Relative  weight  of  terms  (%)  for  density  balance  near  equa- 
torial region:  a and  b for  the  salinity  equation  for  layers 
at  10-m  and  1000-m  depths;  c and  d,  same  as  a and  b except  for 
temperature  equation. 

Same  as  Figure  8 except  for  the  mid-latitude  open  ocean. 
Contour  plot  of  mass  transport  stream  function  in  10^  D 
err? /sec. 

Current  vector  plots:  (a)  at  10-m  depth;  (b)  at  250-m  depth; 
(c)  at  1500-m  depth. 
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Figure  12.  Current  streamline  plots:  (a)  at  10-m  depth;  (b)  at  250-m 
depth;  (c)  at  1500-m  depth. 

Figure  13.  Surface  mean  current  pattern  in  the  North  Pacific  Ocean  (from 
U.S.  Department  of  Commerce,  1961) : (a)  summer;  (b)  winter. 

Figure  Ik.  Vertical  velocity  contours  at  different  levels  in  the  North 
Pacific  Ocean:  (a)  at  20-m  depth;  (b)  at  300-m  depth;  (c) 
at  2200-m  depth. 

Figure  15.  Temperature  distribution  at  different  layers:  (a)  10  m; 

(b)  250  m;  (c)  3000  m. 

Figure  l6.  Observational  surface  temperature  distribution  in  the  North 
Pacific  Ocean  (from  Masuzawa,  1972):  (a)  summer;  (b)  winter. 

Figure  17-  Temperature  profiles  in  the  east-west  vertical  cross  sections: 

(a)  5°N  latitude;  (b)  15°N  latitude;  (c)  25°N  latitude;  (d) 

35 °N  latitude. 

Figure  18.  Temperature  profiles  in  the  cross  sections  from  atlas  or  obser- 
vational data:  (a)  0°-10°N  latitude;  (b)  10°-20°N  latitude; 

(c)  20°-30*N  latitude  (all  from  Moromtsev,  1963) ; (d)  35°N 
latitude  (from  Masuzawa,  1972). 

Figure  19.  Temperature  profiles  in  the  north-south  vertical  cross  sec- 
tions: (a)  105 °W  longitude;  (b)  135 °W  longitude;  (c)  175°W 
longitude;  (d)  155 °E  longitude. 

Figure  20.  Temperature  profiles  in  the  north-south  cross  sections  from 
atlas  and  observational  data:  (a)  100°-110°W  longitude;  (b) 


130°-lU0°W  longitude;  (c)  170o-l80°W  longitude  (all  from 
Moromtsev,  1963) ; (d)  155°E  longitude  (from  Masuzawa,  1972). 
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Figure  21. 


Figure  22. 


Figure  25. 


Figure  2k. 


Figure  25. 


Figure  26. 


Figure  27. 


Figure  28. 


Salinity  distributions  of  layers  at  depths  of:  (a)  10  m; 

(b)  250  in;  (c)  3000  m. 

Average  sea  surface  salinity  in  North  Pacific  Ocean  (from 
Reid,  1969):  (a)  summer;  (b)  winter. 

Salinity  profiles  in  north-south  vertical  cross  sections  at: 

(a)  105 °W  longitude;  (b)  l65°E  longitude;  (c)  l80°W  longitude. 
Observational  salinity  profiles  in  the  north-south  vertical 
cross  sections:  (a)  100°-110°W  longitude;  (b)  l60°-170°W 
longitude  (all  from  Moromtsev,  1963);  (c)  l80°W  longitude  (from 
Moromtsev,  1963). 

Pressure  distributions:  (a)  at  surface;  (b)  at  300-m  depth; 

2 2 

(c)  at  2250-m  depth.  Contour  interval  is  0.6  m /sec  . 

Annual  mean  dynamic  height  in  dynamic  meters  with  respect  to 
1000-m  reference  level  (from  Wyrtki,  197I+). 

Energetic  diagrams  in  the  NORPAX  model.  Units  for  energy  are 

3 

ergs/cm  , and  for  rate  of  energy  transformation, 
ergs/cm^/day. 

Zonally  averaged  meridional  heat  transport  in  the  model. 

H,p  indicates  total  northward  heat  transport;  H , heat  trans- 
port due  to  total  current;  H , , heat  transport  due  to  baro- 
clinic  current;  H-,  heat  transport  due  to  barotropic  current; 
H^,  heat  transport  due  to  diffusion  processes. 


Figure  29. 


Total  meridional  heat  transport  in  the  North  Pacific  Ocean: 

(a)  from  the  ocean  model  (NPO);  (b)  from  Vonder-Haar  and  Oort 
(1973),  (VO);  (c)  from  Budyko  (1963)  and  Sellers  (1965),  (BS). 
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Ocean  energy  transport  estimated  by  Vonder  Haar  and  Oort  (1973) 


